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Abstract 

Q Data assimilation method consists in combining all available pieces 

of information about a system to obtain optimal estimates of initial 
£^ states. The different sources of information are weighted according 

to their accuracy by the means of error covariance matrices. Our 

purpose here is to evaluate the efficiency of variational data assim- 

^__i ilation for the xenon induced oscillations forecasts in nuclear cores. 

In this paper we focus on the comparison between 3DVAR schemes 
)rZ with optimised background error covariance matrix B and a 4DVAR 

00 scheme. Tests were made in twin experiments using a simulation code 

^>0 which implements a mono- dimensional coupled model of xenon dy- 

namics, thermal, and thermal-hydraulic processes. We enlighten the 
very good efficiency of the 4DVAR scheme as well as good results with 
the 3DVAR one using a careful multivariate modelling of B . 
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1 INTRODUCTION 

In this article we aim to evaluate the efficiency of variational data assimilation 
methods known as 3DVAR and 4DVAR methods in the forecasting of Xenon- 
135 oscillations. 

Xenon-135 is known to be at the origin of axial power oscillations of about 
one day period in pressurised water reactors (PWRs)[lJ. These oscillations 
do not change the overall power produced by the nuclear plant but they are 
undesirable from a safety point of view. As soon as oscillations are detected, 
they are damped using appropriate control rod movements inside the core. 
Detection as well as prediction of xenon induced oscillations are an important 
part in the operation of a nuclear power plant. 

No direct measurement of the concentration of xenon in the reactor core 
is available, and the simulation of the nonlinear xenon dynamics still rep- 
resents a challenge for real time applications such as system monitoring. 
Several models have been proposed for the real time estimation of xenon 
concentration. They include flux and iodine-135 dynamics modelling. Some 
of them require an estimation of parameters such as presented in [2] . This 
approach is a first step in improving the state estimation but it does not 
take into account errors in the measurements used to adjust the model pa- 
rameters: therefore bad measurements can affect adversely the quality of the 
computed state. In addition, it does not allow to correct initial conditions 
of the coupled flux-iodine-xenon dynamical system. Most of these models 
assume equilibrium concentrations, though the initial distributions of iodine 
and xenon have a significant impact on the power transient. 

Song and Cho [3] determined an analytic initialisation of iodine and xenon 
of an out-of-equilibrium state which consists in adding a corrective term with 
a sinus shape to the ID equilibrium concentrations. The amplitude of the 
sinus is fitted with axial offset power measurements. These measurements 
are considered to be perfect, in the same way as the xenon dynamics model is 
considered exact. This approach is based on analytical developments which 
limit the pattern of the added correction and still does not take into account 
the errors in the measurements. 

Here we aim to improve the Xenon-135 concentration forecast by finding 
alternative and more accurate solutions using variational data assimilation 
techniques. Such techniques find their root in earth science and are used daily 
in weather forecast. Data assimilation is nowadays more and more used in the 
nuclear science community as well as for the improvement of the nuclear core 
activity determination [U E] as for the nuclear accident model parameter 
determination J6J [7] . Thus it is challenging to apply such a technique to 
xenon oscillation forecast even if other trials have been done through genetic 



algorithm [8] or optimal control [9]. 

Even if the problem we deal with can be seen as very well adapted to the 
Kalman method as it has been shown in [TU], we choose to use the variational 
one instead. This choice has been done as on long term plan we expect to 
use some high resolution models, hopefully together with their adjoint. Thus 
it is important to already test variational techniques that are the only one 
available when the control space is too large. 

Two data assimilation methods, the 3DVAR and the 4DVAR, will be 
used here. The main difference between those methods is that the 4DVAR 
takes into account the dynamic of the process. This is a tremendous im- 
provement as it has been proven in weather forecast. But such a method 
needs the adjoint that is costly to develop in an industrial context. To over- 
come this difficulty a specific model (CIREP1D) compatible with automatic 
differentiation has been developed. 

In a first part (section [2]) we describe the model, CIREP1D, a mono- 
dimensional xenon dynamics model which includes neutron and thermal- 
hydraulic processes. Then we give a brief overview of variational data assim- 
ilation methods in Section [3} The setting of the data assimilation method 
is presented in Section |4} In Section [5] we propose three background error 
covariance matrix modellings (two univariate modellings and a last multivari- 
ate modelling). Finally in Section [61 we compare the quality of the 3DVAR 
estimates to a 4DVAR estimate. 

2 THE CIREP1D MODEL 

Since 3D operational industrial codes are time consuming, we use a monodi- 
mensional axial xenon/iodine dynamics model coupled with a monodimen- 
sional neutronic/thermal/thermal-hydraulic model, named CIREP1D. It sim- 
ulates axial xenon dynamics according to the overall power and the control 
rod insertion records in a given time window. CIREP1D takes a few seconds 
to simulate a xenon oscillation of a one-week time range but contrary to 
simpler models, it gives access to quantities measured in core: axial power, 
axial xenon, axial iodine, axial flux and boron concentration. The agreement 
about axial xenon dynamics between ID and 3D models is good and has 
been studied in detail in [TT] . 

Globally speaking, CIREP1D solves a nonlinear system of ordinary dif- 
ferential equations given by an operator Q: 

d{Cx *f I \ z ) t) = g{c Xe ,c I ){z 1 t) ) (i) 



by using an implicit Euler scheme. Cx e {z,t) and Cj(z,t) represent the axial 
concentration of xenon and iodine respectively. Each time step requires a 
critical boron concentration computation corresponding to the assumption 
that neutron and thermal-hydraulic effects may be treated with stationary 
coupled equations. The xenon dynamics can be initialised by either a given 
xenon and iodine concentrations or by equilibrium concentrations. Hereafter 
we give CIREP ID equations in detail. 

Iodine and xenon equations : Iodine and xenon balance equations are 
differential equations on the variables z, t: 

^ = 7/ E / *-A J C 7> 

(2) 

^p = TX e S/$ + AjCj - (\ Xe + <TXe $) C Xe - 

Sj is the fission cross section of the fuel and a Xe the absorption cross sec- 
tion of xenon- 135. These variables depend on z and t. The z-coordinate is 
measured from the bottom of the ID reactor. 7/ and 7x e are the fractional 
fission yield of iodine and xenon. Finally, A/ and Xx e are decay constants of 
iodine and xenon. 

Neutronic model: The neutronic flux $ = ($i,$2) is identified by solv- 
ing two-group diffusion equations. We assume that the time step of flux 
simulation (a few seconds) is shorter than the xenon oscillations. As a con- 
sequence, at each time step for the resolution of xenon equation, the flux can 
be computed using the stationary diffusion equations: 

-8 z D 1 d z ^ 1 + [S ol + £ r ] $! = jj^E/fc, 

-d z D 2 d z * 2 + Pa2 + D 2 ] $ 2 - S r $x = 0, (3) 

with vT,f& = vi¥,fi<&i + U2^f2^2- 

where $1 et $2 are groupwise neutron axial flux distribution, S r is the scat- 
tering cross section and S a9 , D g and f 9 S/ 9 are the absorption cross section, 
the diffusion coefficient, and the neutron emitted in fission cross section. All 
these variables depend on the z variable. The system is closed by using 
albedo boudnary condition. The balance is obtained by looking for boron 
concentration such that the eigenvalue k is equal to one (this is critical boron 
concentration computation). The boron influence does not appear explicitly 
in the previous equations but is linked to cross section values through the 
feedback model. 



Feedback model: Thermal effects are not lumped in a power feedback 
parameter as done in some other cases[2j. The developed feedback model is a 
linear interpolation model relying on assumption that the cross sections de- 
pend on six quantities: fuel irradiation, xenon concentration Cxe, boron con- 
centration Cb, moderator density p mo d, moderator temperature T mo d and fuel 
temperature Tf. Therefore, CIREP1D includes a thermal/thermal-hydraulic 
model, as described below. 

Thermal-hydraulic model: Since the speed of the water flowing upwards 
through the reactor is high, we can assume that the thermal- hydraulic prob- 
lem is an axial monodimensional problem for the slow transients which are 
common in the normal operational mode. The moderator temperature T mo d 
is then described by the following equation: 

Q d z T mod (z, t) = - [P l / n (z, t) + P l Z d {z, t)} , (4) 

PmodPraod 

where Q, c mo d and p mo d respectively represent the volume flow rate, the 
moderator specific heat capacity and the moderator density. Lineic power 
Pf and Pmod released in both fuel and moderator are computed from the 
known two-group flux $. 

Thermal fuel model: Contrary to the thermal-hydraulic model, we em- 
ploy a radial model for the thermal fuel model. Thus, a radial description of 
the fuel is required. We neglect the axial conduction in fuel pin and assume 
rotational symmetry of the problem. Under these assumptions, the thermal 
problem can be described by a monodimensional model in the radial variable 
r: 

- h f d r T f (r,z) - \ f d 2 r T f (r,z) = P f {z)/A. (5) 

The variable A/ represents the fuel thermal conductivity and A corresponds 
to the pin section. This equation is coupled with the neutron equation 
through the lineic power Pf and to the thermal-hydraulic problem through 
the boundary condition expressed on edge Y: 

Mr e T, X f d r Tf(r, z) = h tot (z) [T f (r, z) - T mod {z)} , (6) 

where h to t is the thermal exchange parameter. The thermal and thermal- 
hydraulic parameters pmod-, c mo d, A/ and h to t depend on moderator and fuel 
temperatures. Therefore, the coupled thermal/thermal-hydraulic problem is 
nonlinear. 



Xenon transient simulation As an example of CIREP1D simulation, we 
present results from a computation with the following characteristics: time 
range of 200 hours, load following during 30 minutes and, then until the end 
of simulation, no further rod movement and no power change. The simulated 
core is in the middle of a burnup cycle and then is moderately irradiated. 
Figure [T] shows a xenon oscillation which disappears without any external 
intervention after 100 hours. 

3 A BRIEF OVERVIEW OF DATA ASSIM- 
ILATION TECHNIQUES 

Introducing data assimilation method in operational simulation has been an 
important step to improve forecasts as for example weather in meteorology. 
The aim is to provide a satisfactory estimation of the unknown true state 
of a dynamical system by combining all pieces of information about the 
system. This information, obtained from measurements (called observations) 
and simulation, is weighted according to its reliability expressed in terms 
of error covariance matrices. In practice, the model gives a simulated state 
called the background state. The purpose of data assimilation is to determine 
a state, called the analysis state, which is closer to the true state than the 
one described solely by the observations or the model. Thus, the analysis 
state can be used to compute a forecast. 

3.1 Concepts and definitions 

We now introduce some concepts and definitions. A discrete model for the 
evolution of physical system from time ti to time i i+1 is described by: 

X(t i+1 ) = M i+ i,i(X(ti)), 

where X and Ai are respectively the model's state vector and its corre- 
sponding dynamics operator. The dynamics M. of the model evolution is 
commonly nonlinear. We note respectively Mjj and Mf the linear tangent 
and the adjoint operators with respect to the vector Xj associated with the 
dynamics model M. between tj and U. The state vector X is usually obtained 
by discretisation of physical fields on a grid. Its dimension is denoted by n. 
The aim is to evaluate the best estimate of the unknown true state, denoted 
X* which is defined by the best possible representation of reality as a state 
vector at an initial time to- The best estimate that we are looking for in the 
data assimilation process is called analysis and is denoted by X a . 



The information about the system that can be used to produce the anal- 
ysis is listed below: 

measurements in the core gathered into an observation vector Y°. Its 
dimension is p. 

the observation operator. The key to data analysis is to take advantages 
of the discrepancy between observations and state vector. Usually ob- 
servation vector and state vector are not defined in the same space. 
They can be compared through the use of a function from model state 
space to observation space called observation operator and denoted %. 
The operator % can be nonlinear. 

an a priori estimate of the true state before the analysis is carried out. 
This estimate is called background state and is denoted X 6 . In most 
cases, the analysis problem is under-determined because observations 
are sparse and only indirectly related to the model variables. The use 
of this background information helps to make it a well-posed problem 
and to introduce some physical knowledge. Usually, this background 
state is generated from the output of a previous analysis. 

uncertainties in the previous data. Background and observation errors are 
defined by: 

e 6 = X 6 -X* and e\ = Y° - K(M,o(X*)). 

The covariance matrices of these errors are denoted by B and R respec- 
tively. Error modelling is a difficult task, mostly because true state X* 
is unknown and the knowledge of the error covariances is approxima- 
tive. But it is a very important step which influences on the quality of 
the analysis. The basic modelling consists in setting up diagonal ma- 
trices where the diagonal elements correspond to the variances of the 
errors on the background or observation vector. When different fields 
are involved in the state or observation vector, it is then possible to 
choose between an univariate or multivariate modelling. In the former, 
the errors between the different fields are assumed to be uncorrelated 
whereas in the multivariate modelling, the errors are assumed to be 
correlated. Usually errors between the different kinds of observation 
are assumed to be independent and then the covariance matrix R is 
diagonal. It is more common to assume correlations for the background 
part but the evaluation of these ones is also difficult. 

In the same way, the analysis state X a is associated to an analysis 
error defined by: e a = X a — X*. And its covariance matrix denoted by 



A is estimated during the assimilation process or as a postprocessing 
procedure. 

Basically, two families of data assimilation methods exist: stochastic and 
variational methods. The most famous stochastic method is probably the 
Kalman filter. This method is still considered as a reference but its appli- 
cation for real data assimilation problems is limited to problems of small to 
medium size due to its huge computational cost involved in matrix computa- 
tion. Several variants of this method have been developed either to reduce its 
computational cost or to remove the assumption on linearity of the used op- 
erators. Variational methods are based on the minimisation of a cost function 
|12j . These methods, 3DVAR and 4DVAR, that can be adapted to nonlinear 
cases and problems of large size, are mainly used in operational meteorology 
and oceanography since the 1990s. Each variational method is equivalent to 
a stochastic filter method under linear assumptions. 

3.2 Variational methods 

We now give here some elements on variational methods. The 4DVAR cost 
function measures the weighted sum of the square of distances J h to back- 
ground state X fe and J° to the observations Y° over a time interval [to 5 ^n] : 

J iDVAR (X(t )) = J b (K(t )) + J: DVAR (X(t )) (7) 

with 

J b (X) = X - [X-X^B- 1 [X-X 6 ] 

n 

JZdvarW = 5 E t Y ° - U(M i>0 (X))] T Il;i [Y° - H(M i>0 (X))] 



i=0 



where weight matrices B _1 and R^ 1 are the inverse of the background and 
observation error covariance matrices at time t{. Minimisation of ([7]) is done 
with respect to initial state X(t )- In practice, the starting point of the 
minimisation algorithm is taken equal to the background X 6 . Evaluations of 
gradient of Jadvar- 

VJadvar (X(t )) = VJ b (X(t )) + VJ 4 ° DVAR (X(t )) 

with 

VJ b (X) = B- 1 [X - X 6 ] 

VJbvAR (X) = - £?=o M^R- 1 [Y? - M(M,o(X))] , 

are required by most minimisation methods which implies that the adjoint 
operator Mf and H T have to be evaluated. 3DVAR method is a cheaper 



alternative to 4DVAR because it does not require the evaluation of the model 
evolution and its adjoint. The 3DVAR cost function is very close to the 
4DVAR one, except for the time sum that disappears: 



J3DVAR(Mt )) = J b mo)) + J£ D VARfr(to)) 
1 

2 



with Ji DVAR (X) =r \ [Y° - H(X)] T K 1 [Y° - U(X)] . 



In a variational assimilation process, the error covariances of the analysis 
can be deduced from the Hessian of the cost function Jzdvar or Jadv ar\^-?>\'- 

n 

A 4DVAR = B ~ L + J2^ HMi '°^2 DVAR R ' l ( HM ifi)\Xt DVAR - 
i=0 

3.3 Twin experiment frame 

To validate assimilation schemes independently of the model, one performs 
twin experiments. In twin experiments, the initial true state X*(0) is choosen 
and the true trajectory for any time is known. It is a simulated state usually 
obtained by the model used for assimilation. Twin experiments also offer 
the opportunity to compare the analysis to the true state. The true state 
can be used to build background state, for example by adding a noise to X*. 
It can also be used to build synthetic observations by applying observation 
operator % to X 4 and noise afterwards. 

4 COMPONENTS OF THE ASSIMILATION 

SYSTEM 

We develop two variational schemes 3DVAR and 4DVAR, in order to improve 
the xenon and iodine concentration estimation in core, in twin experiments 
set-up. In this section we describe all the components of the assimilation 
system except the B matrix modelling that will be discussed in detail in 
section [5j 

Model 

The evolution model corresponds to the xenon dynamics model imple- 
mented in CIREP1D. This model is based on the resolution of the xenon and 
iodine mono-dimensional time equations. Each iteration time step requires a 
critical boron concentration computation which includes successive station- 
ary neutron/thermal/thermal-hydraulic computations. For such a computa- 
tion, CIREP1D inputs are: 



• initial and final times to and t n of transient, 

• initial xenon and iodine concentrations at time to, 

• transient data: overall power and control rod position variations over 
the time interval [to;t n ]. 

State vector 

The state vector X corresponds to xenon and iodine axial concentrations 
discretised on the 30 nodes of the ID axial spatial mesh used in CIREP1D. 
The dimension of X is then 60. The analysis problem is to find a correction 
<5X such that X a = X b + 5X is as close as possible to X*. This correction is 
searched in the same space as the state vector one. Thus, the minimisation 
problem has dimension 60. 

The observation operator % 

In the present use of data assimilation the observation operator % is 
given by the model itself. It is non-linear and roughly corresponds to a 
critical boron calculation with CIREP1D. Therefore, it depends on xenon 
but not on iodine as no time evolution are done in such a calculation. Since 
the 3DVAR scheme does not involve any evolution model, it cannot control 
iodine concentration. Another important characteristic of this scheme, is the 
quasi-equivalence in computational cost of evaluation of the model M. and 
observation operator H. 

Observations 

In this framework, observations used further in the analysis process are 
not coming from real core measurements. They come from numerical simu- 
lations with CIREP1D, in a twin experiment framework as described in the 
section [3j The scheme is the following: 

1. We compute xenon dynamics initialised by equilibrium concentrations, 
in a time range of one hour for example. The concentrations obtained 
after this hour are defined as the real state X* at the initial time to of 
the future analysis process. 

2. We do a reference simulation with CIREP1D to make the real state X* 
evolve from t to t n : 

X'fe) = M,o(X(t )). 

3. Observations over time range are obtained by introducing a measure- 
ment noise e on real data: 

10 
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The observation vectors Y° bs at different observation times tj are com- 
posed of 3 different measurements: 6 integrated powers over several cells 
(index P), 1 power axial offset data (index AO) and 1 boron concentration 
data (index (7b). The associated error standard deviations are denoted re- 
spectively by aRp, <Jr ao and <jRc b - The observation vector has dimension 
8. 

Error covariance matrices 

To build the observation error covariance matrix R, we assume that mea- 
surement errors are Gaussian, that they are not correlated in space, and that 
they does not depend on time. In this case Ri = R and R is diagonal and 
given by: 



(8) 



Errors introduced in measurements are set to 10%, 5% and 1% of the 
current value respectively for axial power (aR P ), power axial offset (&r ao ) 
and boron concentration {<jRc b ) measurements. Those values correspond to 
the typical knowledge we got on those measurements considering both their 
intrinsic error and the representativity error. 

Minimisation 

Finally, to solve the non-linear minimisation problem, we use the quasi- 
Newton method LBFGS [Hj. This method requires the computation of the 
gradient of J which is done using the adjoint of the xenon dynamics model. 
In our case, the adjoint is obtained by automatic differentiation of CIREP1D 
using TAPENADE software [T5]. Both J and Vj7 are computed in the 
framework of the PALM assimilation coupler [T6l ITT] . 

5 BACKGROUND ERROR COVARIANCE 
MATRIX MODELLING 

The B matrix is one of the most important point of the data assimilation, 
in particular for the 3DVAR method. For the 4DVAR method, the matrix 
is less crucial since the model itself contributes to spread the information. 
Thus in order to have a reliable comparison between both methods a careful 
study of B is presented here. 
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We recall that the state vector of size 60 is composed by the values of 
the physical fields xenon and iodine at the 30 mesh nodes (finite differences 
discretisation), the mesh node numbering starting from the bottom of the 
core. The vector e b of size 60 represent the background error made on each 
of the 30 mesh nodes, which is assumed Gaussian. The error covariance 
matrix B is defined by: 



B = E 



(e»-E[e b ])(e»-E[e<>])- 



We study three different types of modelling for the background covari- 
ance matrix B: an elementary modelling where B is diagonal, a univariate 
modelling where B is block diagonal and at last a multivariate modelling. 

5.1 Settings used for the modelling 

In what follows, we need to use the evolution model M. and then to set up 
a simulated case. We are working in the twin experiment framework. The 
simulated case is a regular transient where the produced power is close to 
the nominal power and control rods are partially inserted in core. The state 
of the core corresponds to the end of a fuel cycle. The true and background 
trajectories X 4 (t) and X b (t) are computed with non equilibrium initial states 
issued from close but not identical previous calculations. 

We use standard deviations axe and 07 close to 3% for the initial B 
matrix. Those values can evolve respect to the treatment we do on this initial 
matrix. Those values come from comparison between CIREP1D and other 
models. Moreover they are in the typical range of values used in JTSJ SI IS]- 

5.2 Elementary modelling 

As a first step to build the B matrix, we omit correlation in space and 
between species consider the diagonal matrix given by: 



B, 



(<&. 
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5.3 Univariate modelling 

Before developing a multivariate modelling, we propose to take into account 
spatial correlations for xenon and iodine. We are looking for a block diagonal 
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matrix: 

B Xe 



B " l Bj 

where the block diagonals are given by: 

Bxe = Bd,XeT'Bd,Xe, 

B/ = B dI TB dI . 

The matrices Bd,Xe and Bd,i correspond to the sub-matrices extracted 
from Bd. The matrix T is built thanks to the Balgovind correlation |19j 
between two nodes of the spatial mesh Zj and Zj numbered from 1 to 30, 
which reads: 

T(zi,Zj) = (1 + \z,i - Zj\/L)exp(-\zi - Zj\/L) 

where the parameter L corresponds to the correlation scale to set up. 

This modelling assures the definite positivity of B M . The choice of L has 
consequences in: 

• the structure of B u (decay property of the matrix elements away from 
the main diagonal), 

• the conditioning of B u (the larger L is, the worse the conditioning is), 

• and the quality of the analysis. 

In practice, the choice of L = 4 for both species gives satisfactory results 

mi. 

5.4 Multivariate modelling 

We propose to build correlations thanks to the evolution model -Mi,o between 
times to and t{. This method is very close to what is done in Kalman filter. 
However we do no consider each step of evolution. 

If we consider a small initial perturbation e on the state X at time t 0) one 
can write: 

Mi,o(X(*o) + efo)) w Mi,o(X(to)) + M i]0 | X .e(t ), 

where M i0 |x represents the tangent linear of .M^o with respect to the 
vector X. Then, if we set e(^) = Mi,o(X.(t ) + e(t )) - Mi,o(X.(t )) , the last 
relation can be expressed in terms of the error as follows: 
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e(U) « Mi )0 |x.e(t )- 

Then we get an approximation of the covariance matrix at t{ as a function 
of the covariance matrix at t : 

cov(e(ti), e(U)) « M im cov(e(t ), e(£ ))Mj 0|x . (9) 

It is noted that it can be difficult to get the tangent linear operator M for 
an industrial code, since it usually requires to be written at the same time as 
the direct code. Thanks to Eq. (J9J) we can model correlations between xenon 
and iodine by multiplying an univariate matrix B n by Mj |x and M^ , x : 

Bj = M ii0 |xB u M ii0 | X . 

In what follows, we compare the univariate matrices B^, B u to the mul- 
tivariate matrices B 3 , B 12 and B 2 4- 

5.5 Overview of the estimated covariance matrices 

Figure [2] shows the diagonals of the various B matrices. These diagonals 
correspond to the error variance of the background state vector for each 
node of the spatial mesh (x-axis in the Figure [2|. We may notice: 



In the diagonal and univariate modellings (Figures 2a and 2b ), the 
variance is the same for all the spatial discretisation points, though 
it is known that the lower and upper regions of the core are quite 
inaccurately modelled. The iodine error variance is bigger than the 
xenon one but it is due to the fact that iodine concentration in core is 
bigger than the xenon concentration. In terms of relative errors, they 
are similar. 



For multivariate modellings (Figures 2c, 2d and 2e), 3 regions can be 
seen (easier to see for the iodine part than for the xenon part): two 
large regions for half lower and upper parts of the core and a small one 
for the central part of the core. The variances are lower for the median 
part than for the two others (also easier to see for the iodine than for 
the xenon) 

The variances are lower in the multivariate modelling than in the uni- 
variate modelling, that is to say tr(B mu iu) < tr(B u ). But the multi- 
variate modelling takes into account correlations between species: thus 
error statistics are spread all around the elements of the covariance 
matrix. 
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To get an idea of the diffusion of the xenon error to the iodine error in the 
multivariate modelling, Figure [3] gives the absolute values (for the clarity of 
the figure) of the correlations with respect to the nodes of the spatial mesh, 
where the correlations are defined as 



corr(i,j) = B(i,j)/y/B(i,i)B(j,j), 

where the element (i, i) of the matrix B corresponds to the variance of the 
background error on xenon at the node number i (resp. iodine at the node 
number i — 30) if i < 30 (resp. i > 31) and the element (i, j) (or (j,i) since B 
is symmetric) corresponds to the covariance of the background error between 
xenon and iodine at the node i. Then the range of the scale varies from to 
1. 

We may notice: 

• Compared to the univariate modellings for which extra block diagonal 
terms are obviously null, introduction of spatial correlations fills in the 



two diagonal blocks in Figure [3bJ The choice of the correlation scale 
will be discussed later. 



In the multivariate modellings (Figures 3c, 3d and 3e), the correlation 
matrices present an internal structure more or less pronounced with 
respect to the use of an evolution model on a more or less long time 
range. As for the variances, one can see two large regions for the half 
upper and lower parts of the core and a narrower region for the central 
part of the core. 

— Diagonal blocks: they correspond to the spatial correlations for 
a given species. Spatial correlations increase with the length of 
the time range used in the evolution model, first for the iodine 
correlations (e.g. Aiu) and then for the xenon ones. This can be 
explained by the way of production of xenon which is produced by 
radioactive decay of the iodine. For the matrix B24, correlations 
between upper and lower regions of the core are almost as strong 
as the spatial correlations inside these regions. 

— Extra-diagonal blocks: they correspond to the correlations be- 
tween the xenon and the iodine. The same correlation increase is 
noticed with the length of the time range used in the evolution 
model. But these correlations stay below the spatial correlation 
level inside the same species. 

— Central blocks: a region including very few nodes of the middle 
of the core seems to be insensitive to the time range used in the 
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evolution model: the background error for these nodes is slightly 
correlated to the background error of the other nodes. At last, 
one can notice that the correlation scale for these nodes seems to 
be shorter to the scale 4 set up in the correlation scale modelling 
described in Section 15.31 

6 COMPARISON BETWEEN 4DVAR AND 
3DVAR SCHEMES 

We are in the twin experiment framework: the true state X* is a simulated 
state. The simulated case was briefly described in the previous section. We 
always assume that background and measurement errors are Gaussian. Er- 
rors introduced in measurements are set to 10%, 5% and 1% respectively for 
axial power, power axial offset and boron concentration measurements. 

The purpose is to compare the result of the 3DVAR and 4DVAR schemes 
at the assimilation time but also after at a forecast of 10 hour (typical time for 
xenon oscillation). In the case of the 3DVAR assimilation, various modelling 
of the B matrix will be under consideration. The 3DVAR results are com- 
pared to a 4DVAR result whose characteristics are set up according to [20J: 
the window size is set to 6 hours and the observation frequency is set to 2 
hours (then 3 observation sets are used) 

The various B matrices used for the 3DVAR schemes are the following: 

• Bd corresponds to the univariate modelling of B without any spatial 
correlation; 

• B M corresponds to the univariate modelling with spatial correlation; 

• Bj corresponds to the multivariate modelling where the time range used 
in the evolution model equals to i hours (here 3, 12 or 24 hours). 

Results are organized as follows: firstly we are analysing the 4DVAR 
results. Secondly we are comparing ID errors on axial shape of xenon, iodine 
and power. Then we show the time evolution of the mean error. At last we 
discuss the statistics of the analysis given by the analysis matrices. 

6.1 4D VAR results 

First we will look at the results of the 4DVAR scheme that is known to 
be the most efficient for forecasting. With this scheme, the modelling of 
xenon/iodine correlation for B has a weaker influence on the quality of the 
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analysis and then we use the B u matrix. Figure |4] shows relative error^] on 
xenon, iodine at the assimilation time to an d on power at to an d t + lOh for 
two computations: the background computation which gives the background 
state and a 4DVAR computation with the assimilation characteristics given 
before. 

On this Figure [4j one can see that the 4DVAR computation allows to 
reduce errors on xenon and iodine axial shape at least by a factor of 2. For 
the axial power shape the decrease is even more important since the errors are 
reduced by a factor of 4. What follows aims at showing that it is possible to 
approach the 4DVAR results quality with 3DVAR schemes using multivariate 
B matrices. 

6.2 Comparison between 4DVAR and 3DVAR 

Figure [5] gives the ID relative errors respect to the true state on axial shape of 
xenon, iodine and power for the different assimilation schemes. We see that 
the 4DVAR results are better than the 3DVAR ones independently of the 
choice of B. However, looking at 3DVAR results at analysis time, we notice 
that xenon is rather well estimated in all variational schemes except for the 
scheme which uses the elementary modelling of B. The xenon estimation 
is not as good for the bottom half as for the top half, but it is everywhere 
better than the state given by the background computation. This good 
result can be explained by the fact that the xenon level at the assimilation 
time is directly related to the assimilated observations, that is to say to the 
integrated powers. 

On the opposite, iodine is not directly related to the power level but 
through the production of xenon since xenon is essentially produced by the 
radioactive decay of the iodine. Thus the observation operator Ti does not 
depend on iodine. Therefore it is not possible to correct iodine state with a 
3DVAR scheme unless correlations between xenon and iodine are introduced 
in B. As a consequence the 3DVAR analysis error is equal to the background 
one for the computation with the matrices B^ and B u . With a multivariate 
modelling of B it can be seen that the time range used to build the extra- 
species correlations influences the quality of the iodine analysis. It seems 
that the longer the time range is, the better the analysis is. In fact there is 
no convergence towards the 4DVAR analysis quality, and taking a 48 hours 
time range does not allow to improve 3DVAR results. One can assume that 
the optimal time range is related to the time constants of radioactivity decay 



lr The "true" value of the fields xenon, iodine and power is known since we are working 
in the framework of twin experiments. 
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of the xenon and iodine. 

The knowledge of the iodine level is not important for the monitoring 
operator system. We are expecting to improve axial power forecasts with the 
assimilation techniques. If we take a look at the power estimation after 10 
hours, we notice that background and diagonal and univariate 3DVAR anal- 
ysis axial power shapes are very close. Since xenon is essentially produced 
by radioactive decay of iodine, a bad estimation of the initial iodine con- 
centration will affect the xenon concentration estimation later. As expected 
the multivariate modelling leads to a significant improvement in the 3DVAR 
forecast (up to a factor of 2 to 3 on the errors for the use of the matrices B12 
and B24). The multivariate modelling resulting from the use of a time range 
of 3 hours in the evolution model does not seem to be sufficient to really 
improve power forecasts. As a first conclusion, the optimal time range seems 
to be around 12 hours: the 24 hours time range does not improve much the 
forecasts compared to the 12 hours time range but it increases the amount 
of computation for B. 

6.3 Time evolution of the L 2 -norm error 

We would like to confirm what has been shown at the assimilation time to 
and to + 10 hours. Figure [6] presents the time evolution of the L 2 -norm error 
of the three studied fields for the different matrices B used. An oscillation 
which tends to damp can be seen. But it does not change the previously 
drawn conclusion. Two groups can be seen: the first one composed by the 
analysis computed by the 3DVAR schemes using the matrices B^, B u , B 3 ; 
the second one composed by the analysis computed by the 3DVAR schemes 
using the multivariate matrices B i2 , B 2 4 plus the analysis computed by the 
4DVAR scheme based on the matrix B n . As one can see on the Figure [6j the 
latter still represents a "reference" towards one can except to tend with a 
3DVAR scheme based on a multivariate modelling. And the 3DVAR schemes 
where no or very few correlation between xenon and iodine is modelled are 
unable to give good forecasts. We can deduce that the 3DVAR assimilation 
results depending a lot of the B involved . 

6.4 Analysis matrices 

Assimilation techniques offer a posteriori diagnostic on the computed analysis 
by the mean of the analysis matrix. We use this opportunity to go more 
deeply into the study of the respective quality of the various approaches. 

Figures [7] and [8] are to be compared to Figures [2] and [3] that show the 
structure and amplitude of the elements of the background matrix B. Figures 
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7a] and [7b] show that the 3DVAR schemes based on univariate process is not 
able to reduce the variance (represented by the diagonal of A) of the iodine 
error. Though the scheme based on the multivariate matrix A 3 has been 



proved not being good for forecast, Figure [7c] shows that it is much better 
than the previous ones since it allows to reduce by a factor of 2 the error on 
iodine. 

Figure [8] gives absolute values of the correlations with respect to the nodes 
of the spatial meshes. It is to be noted that if the matrix B does not contain 
correlations between the xenon and iodine errors (extra diagonal blocks), 
the analysis matrix A issued from a 3DVAR scheme does not contain any 
correlation between them neither. The 4DVAR scheme is based on the uni- 
variate B u but the corresponding analysis matrix A^ovar shows correlations 
between xenon and iodine. These latter are brought by the 4DVAR scheme 
through the use of the evolution model on a 6 hours time range but they are 
shorter than the ones in the matrices A3, A 12 and A24. 

Figure [8] also presents the spatial correlations intra species (diagonal 
blocks). It can be compared to the spatial correlations shown in Figure 
[3] It can be seen that spatial correlations after the assimilation process are 
a little bit shorter than they were in the matrix B. 

7 CONCLUSIONS 

In this paper, we have shown how variational data assimilation methods can 
be used to improve the accuracy of the prediction of the xenon concentra- 
tion in PWRs. A mono dimensional xenon dynamics code CIREP1D was 
developed for this purpose. 

The investigation done here in twin experiments proves that the 4DVAR 
scheme is a very efficient method to improve the accuracy of the prediction of 
the xenon concentration as well as the axial power shape in PWRs. However 
this method is computationally costly and the development of the adjoint of 
the model is mandatory. 

A computationally cheaper solution is the 3DVAR that can lead to rather 
good result. Nevertheless these latters can only be obtained through a careful 
modelling on the associated background error covariance matrix B. Among 
the various modellings of B studied here the one based on the multivariate 
modelling is the most satisfactory. 

The next stage for the 3DVAR approach could consist in setting up an 
assimilation chain where the B matrix is updated at each analysis step as it 
is successfully done for meteorological operational forecasts. 
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Figure 1: Xenon dynamics simulation on a time range of 200 hours. The 
map of the first figure represents the evolution of the monodimensional xenon 
vector with respect to the time (x-axis). The y-coordinate corresponds to 
the position in core and the vertical coordinate gives the xenon level. The 
second plot gives the xenon axial offset with respect to the time, i.e. the 
power difference between the half top and bottom parts of the core. 
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Figure 2: Diagonal terms of B matrices for univariate and multivariate mod- 
elling, that is to say background error variances with respect to the node of 
the spatial mesh. 
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Figure 3: Absolute values of the components of the correlations related to 
the matrices B for univariate and multivariate modelling with respect to the 
nodes of the spatial mesh. 
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Figure 4: Relative errors on xenon, iodine at to an d on power at to and 
to + lO/i for the background state and the 4DVAR analysis state. 
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Figure 5: Relative errors on xenon, iodine at to and on power at to + lO/i. At 
to 3DVAR iodine analysis and background are mistaken. A to + lOh, 3DVAR 
and background predicted axial power shapes are mistaken. 



27 




12 IS 24 30 36 
(hours) 



Figure 6: Relative L 2 -norm errors on monodimensional xenon and iodine 
fields with respect to time for 3DVAR and 4DVAR schemes based on various 
background matrices. 
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Figure 7: Diagonal terms of the analysis matrices A for 3DVAR and 4DVAR 
schemes based on univariate and multivariate modelling of B. 
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Figure 8: Absolute values of the components of the correlations related to 
the matrices A for 3DVAR and 4DVAR schemes based on univariate and 
multivariate modelling of B. 
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